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We give a general scheme for finite-differencing partial differential equations in flux-conservative 
form to second order, with a stencil that can be arbitrarily tilted with respect to the numerical 
grid, parameterized by a "tilt" vector field j A . This can be used to center the numerical stencil 
on the physical light cone, by setting ~f A = (3 A , where (3 A is the usual shift vector in the 3+1 split 
of spacetime, but other choices of the tilt may also be useful. We apply this "causal differencing" 
algorithm to the Bona-Masso equations, a hyperbolic and flux-conservative form of the Einstein 
equations, and demonstrate long term stable causally correct evolutions of single black hole systems 
in spherical symmetry. 



I. INTRODUCTION 

The apparent horizon boundary condition (AHBC) ap- 
pears to be one of the fundamental techniques required 
for evolving black hole spacetimes using numerical tech- 
niques. In this paper, we present work on an AHBC in 
the context of the recently formulated Bona-Masso (BM) 
hyperbolic system for the Einstein Evolution equations, 
however in doing so, we present a technique which is gen- 
erally causally correct for any first order flux-conservative 
set of PDEs. 

The idea of the AHBC was credited to Unruh by 
Thornburg H. The fundamental idea is that, rather than 
avoiding a singularity by taking slices which delay the in- 
fall of observers inside the horizon (eventually requiring 
an infinite force), one could take regular slices everywhere 
outside the horizon and place some stable boundary con- 
dition inside the apparent horizon, excising a portion of 
the numerical grid. Mathematically this is consistent be- 
cause the apparent horizon is known to be inside the 
event horizon, and the interior of the event horizon is, by 
definition, causally disconnected from its exterior. 

The AHBC in numerical relativity was first shown to 
work by Seidel and Suen in Ref . j| . The AHBC proposal 
by Seidel and Suen has two components. The first is to 
choose a shift condition which, after some evolution, locks 
the coordinates by freezing the position of the horizon 
and keeping radial distances between points constant. In 
spherical symmetry, this uniquely determines the shift 
everywhere. Additionally, to handle large shift terms, 
Seidel and Suen propose re-writing the finite difference 
representation of the ADM equations to obey the causal 
structure of the spacetime. 

This causal differencing is an important aspect of much 
AHBC work to date and is generally credited to Sei- 
del and Suen ("causal differencing") or Alcubierre and 
Schutz ("causal reconnection" ) j3|. The idea proposed 
by Seidel and Suen is as follows. In the presence of a 
shift, the Einstein equations have additional terms in the 
evolution equations (due to the action of Cp on jij and 
Kij). However, there is another coordinate system in 



which the shift is zero. Finding the coordinate system, 
in general, involves integrating a transformation function 
in time. The Seidel and Suen prescription is to finite 
difference in the transformed zero shift co-ordinates and 
then re-transform the finite difference representation to 
coordinates with a shift. 

Using these two techniques, Seidel and Suen pro- 
ceed to demonstrate that they can accurately evolve 
Schwarzschild black holes and Schwarzschild black holes 
with infalling scalar fields for long periods of time. 

The followup to Seidel and Suen, a paper by Anni- 
nos, Daues, Masso, Seidel and Suen, j| gave details of 
the Seidel and Suen causal differencing scheme and pre- 
sented several shift choices. These shifts allow for co- 
ordinate regularity in the entire spacetime and provide 
horizon locking. Moreover, several of these shifts are 
extensible to full three dimensional cases, most notably 
the minimal distortion shift Using causal differenc- 
ing and horizon locking shift conditions, Anninos et al. 
are able to evolve Schwarzschild black holes in spherical 
symmetry for 1000M with very small errors in the mea- 
sure of the mass of the horizon, compared to 100% mass 
errors present in simulations without an AHBC around 
t = 100M. 

The first application of an AHBC to a hyperbolic 
scheme was that of Scheel et al. Q . This scheme used a 
hyperbolic formulation due to York on a Schwarzschild 
black hole. The essence of the causal difference scheme 
was to decompose the time derivative into time evolution 
along the normal direction and spatial transport due to 
the shift. That is, they would evolve along the normal 
direction to some point no longer on their numerical grid, 
and then re-construct their numerical grid by interpola- 
tion. The Scheel et al. approach is roughly equivalent to 
our advective (non flux-conservative) causal differencer 
with interpolation after the step, as discussed below. The 
results of Scheel et al. were initially disappointing, as an 
instability arose on a short (10— 100M) time scale. Very 
recent work J7j] has removed this instability by adding 
constraints to the evolution equations (in a manner spe- 
cific to spherical symmetry), and run times exceeding 
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lOOOOAf have been achieved. 

Anninos et al. extended their methods by importing 
a one dimensional shift into a three dimensional code in 
||. They found a stable horizon mass for a moderate 
time evolution. This work was extended in the thesis of 
Greg Daues ||, who applied various shift and excision 
conditions to three dimensional systems evolved in the 
ADM formulation, successfully evolving a Schwarzschild 
hole for 100M using live gauge conditions in three dimen- 
sions. 

Another demonstration that some aspects of an AHBC 
could work in three-dimensional numerical relativity was 
given by Briigmann [ pl)[ . Brugmann does not use a shift 
to freeze the horizon, and only demonstrates his method 
in geodesically sliced spacetimes. That is, the horizon 
keeps swallowing grid points as the simulation continues. 
Nonetheless, he uses an irregular inner boundary on a 
(semi-adaptive) Cartesian grid and demonstrates a stable 
evolution for several times longer than the irM infall time 
of the throat in a geodesically sliced spacetime. The inner 
boundary is simply a few zones inside the location of 
the (moving) surface r = 2M, rather than a numerically 
located horizon. The inner boundary is generated by 
polynomial extrapolation. 

Recently, the Binary Black Hole Grand Challenge has 
presented several convincing results using an apparent 
horizon boundary condition with an ADM system using 
a causal differencing scheme similar to that of Scheel et al. 
extended to three dimensions. Using boosted Kerr-Schild 
slices of Schwarzschild with exact gauge conditions, the 
Grand Challenge was able to transport a black hole across 
a grid . The Grand Challenge also reports the ability 
to hold a static black hole static for close to 100M in 
Eddington-Finkclstein slicings in three dimensions. 

The work presented below is an attempt to apply sim- 
ilar techniques to those used by Daues and by the Grand 
Challenge to the BM evolution system, while attempting 
to exploit the first-order, flux-conservative form of the 
BM evolution equations. In the first half of the paper, 
we lay a groundwork for future extensions of the BM sys- 
tem to general three-dimensional BH spacetimes, which 
are far beyond the scope of this paper, by discussing 
and analyzing various fully three-dimensional techniques 
for implementing a general AHBC. In the second half, 
we apply the general framework to the BM system ap- 
plied to spherically symmetric vacuum spacetimes, that 
is, to the Schwarzschild black hole. We evolve initial data 
taken from three different time-independent slicings of 
Schwarzschild, all of which are regular at the horizon. 
We use a lapse and shift imported from the exact solu- 
tion on all of them, and an exact shift together with a 
live harmonic slicing condition on one of them. In each 
of these situations we test one traditional and four causal 
finite differencing schemes. 



II. CAUSAL DIFFERENCING OF 
FLUX-CONSERVATIVE EQUATIONS 

In this section, we develop causal differencing for an ar- 
bitrary system of flux-conservative equations. The equa- 
tion we want to finite difference is 



du dF A {u) 
dt dx A 



S(u), 



(1) 



where x A are the spatial coordinates, and u is a vector 
of unknowns. A necessary condition for a stable finite 
differencing scheme is that all the characteristics of the 
system lie inside the numerical domain of dependence 
("the stencil"). Depending on the formulation of the 
Einstein equations, there may be mathematical charac- 
teristics other than, and in fact outside, the physical light 
cone. In this paper, we consider only situations where the 
mathematical characteristics lie inside the physical light 
cone. In Fig. |l| we show the relationship between the 
numerical and physical light cone, and show situations 
in which stable or unstable situations result including a 
causal differencing case, where we adjust the numerical 
stencil to follow the physical light cone. The character- 
istics are not immediately apparent from the form (|l|) of 
the equation, but if we are dealing with a hyperbolic for- 
mulation of the Einstein equations, perhaps coupled to 
matter, we know that the (physical) characteristics are 
centered on the center of the light cones of the spacetime 
metric 

ds 2 = -a 2 dt 2 + g A B{dx A + (3 A dt){dx B + (3 B dt). (2) 

We now introduce an auxiliary coordinate system (t,x a ) 
that is related to (t, x A ) by 



where x a (t,x A ) obeys the differential equation 



dx a A dx a 



dt 



7 



dx A 



(3) 



(4) 



The vector field j A distorts the x coordinate system rel- 
ative to the x coordinate system. We define the short- 
hands 



„ di a dx a 
m = M a A - 



dt ' 



dx A ' 



Note that because t = t, Q is also the matrix inverse of 
M. 
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A. Transforming the differential equation 



t+dt 



t+dt 





No Shift, Non Causal 
Stable 



Shift, Non Causal 
Stable 



t+dt 



t+dt 





Shift, Non Causal 
Unstable 



Shift, Causal 
Stable 



FIG. 1. We show the stability criterion for light cones and 
numerical stencils. The light gray areas show the physical past 
light cone of a given point, and the dark gray areas show the 
numerical stencil. In the top two figures, where the physical 
light cone is inside the numerical light cone, we will get a 
stable evolution. In the non-causal large shift picture, the 
evolution will not be stable, but tilting the numerical light 
cone (or causal differencing) will result in a stable evolution. 



Assuming for a moment that we choose j A — (3 A , 
spacetime metric in the new coordinate system is 



ds 2 = -a 2 dP 



g a b 



dx a dx 



7,b 



the 



(0) 



This has no shift, and therefore the light cone is symmet- 
ric around d/dt. The 3-metric, lapse and shift in the two 
coordinate systems are related by 



g AB =M a A M b B g ab , 



a, 



(3 A = Q J 



(7) 



Our motivation for introducing the auxiliary x coordi- 
nates is causal differencing. Here we take causal differ- 
encing to mean using a stencil that is symmetrically cen- 
tered on the light cone. This is equivalent to its being 
centered on the vector field d/dt in the new coordinate 
system. We obtain a "tilted" stencil in the x coordinates 
by first transforming the differential equation to the x co- 
ordinates, using a finite differencing scheme to advance 
in x coordinates, and transforming the result back to the 
x coordinates. Although 7^ = f3 A is the choice that mo- 
tivates our scheme, we should keep in mind the point of 
view that 7 simply parameterizes a family of tilted sten- 
cils for the differential equation ([!]). No reference needs 
to be made to either the spacetime metric or the fact that 
([!]) is the Einstein equations. In this respect our scheme 
differs from the causal differencing schemes suggested by 
Seidel and Suen |^| and Scheel et al. Q , which explicitly 
use the fact that the Einstein equations simplify in the x 
coordinates. In the following sections we consider again 
a generic tilt vector field j A . 



The partial derivatives in the two coordinate systems 
are related by 



(8) 
(9) 
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Transforming (Q) to (t, x a ) we obtain 
„ du 



du 
~dt 



dx a 



QpA 



(10) 



We note that this equation is not in flux-conservative 
form. We will refer to ( |Io|) as the advective or non flux- 
conservative form of our evolution equation. We can put 
( |To| ) in flux-conservative form by moving m a and M a a 
into the fluxes, giving additional source terms. Doing 
this, we obtain 

% + ("»"« + MU aF A ) = S(u) + A u + A A F A , 



(11) 



where the flux correction coefficients m a and M a A have 
already been defined, and the source correction terms are 



d d 
A A = -7p—M a a = t—t lndet M, 

dx a dx A 



X 



d 
dx a 



d_ 

dt 



In det M. 



(12) 
(13) 



The second equalities have been written down to clarify 
the function of the source correction terms, namely to 
account for a change in the volume element detM. It is 
useful to note the identity 



A 



d~f A 
dx A 



■-y A A A . 



(14) 



We emphasize that when applied to the Einstein Equa- 
tions, S(u) and F A are the BM sources and fluxes includ- 
ing the shift terms. Even in the case of r y A = (3 A , we do 
not cancel these terms from the sources and fluxes, but 
allow the cancellation to take place numerically (both in 
the advective and flux-conservative schemes). Below we 
will explore only the j A = (3 A case numerically, but will 
find interesting theoretical advantages in the j A = t[3 a 
case (where r is some constant). 

The advective form of the equation, (|l0|), has several 
potential numerical advantages over the flux-conservative 
form, ([n]). Most notably, the equation d t u — for u = 
constant requires cancellations between flux derivatives 
of the shift and divergences of the shift in the sources in 
(|il"|), while no such cancellations are needed in (|l0|). 

We make the two coordinate systems agree at t = t = 

0: 
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$ a (x A ,Q) = 6 a AX / 



(15) 



In the remainder of this work we also assume that the 
tilt is "frozen" in the x coordinates: 



d_ 

at 



7 



0. 



(16) 



In practice, the tilt will be constant only throughout one 
time step, as described in section II C. For convergence 



tests, we shall choose it to be constant throughout one 
time step on the coarsest grid, and accordingly several 
time steps on the finer grids. After each time step we 
discard the coordinate system (t, x a ) and start again from 
scratch. We identify the two coordinate systems either 
at the end or at the beginning of the time step. The 
implication of this choice on boundary conditions will be 
discussed below. We will label the time when the grids 
coincide as t = 0, so time steps will go from t — —At to 
t = or from t = to t = At. 



-7 C (7^),c, 

B A . . _ „,B , 



E (1)B 
A(i)A 

A(2)A = 7,A A (1)B - l"h-(l)A, 



B ■ 



(23) 
(24) 
(25) 



Note that all the Taylor coefficients on the right-hand 
sides are evaluated at t = t = 0, where the x and x grids 
coincide. As furthermore j A is independent of time on 
points on the x grid, we can think of the right-hand sides 
as simply evaluated on the x grid (at whatever time). 
For the same reason, one obtains these quantities on half- 
points of the x grid (at any time) simply by averaging the 
Taylor coefficients between neighboring points on the x 
grid: the result is already accurate to quadratic order. 
Note also that A and m a are exactly constant along x 
lines. 



C. Treatment of the lapse and shift 



B. Calculating the source and flux correction terms 

We need to know either the x A traced back to t = — At 
along lines of constant x or x a traced forwards along lines 
of constant x A , in order to interpolate the data at the 
beginning of the time step onto the x grid or reconstruct 
the x grid after the time step. Rather than by solving 
the differential equation (|J) by finite differencing, we do 
this by a Taylor series expansion in t: 



x A (x, t) = S A a x a -j A t+ -7 B 7j73 t 2 + 0(t 3 ), (17) 



x a (x,t) = S 11 , 



x A + 1 A t + ^ B j A t 2 + 0(t 3 ) 



(18) 



In order to calculate the corrected fluxes and sources in 
the differenced version of ( |Tl] ) , we also need to know m a , 
M a A, A and A^ on the i a grid, at some of the times t = 
±At, t — ±At/2 and t = (exactly which times depend 
on our numerical integration scheme; Lax-Wendroff will 
require the half-times, MacCormack will require the full 
times, and so forth). [For ( [To|) we need m a and M a A .] 
There are different ways of evaluating these quantities. 
We have chosen to expand all auxiliary fields in a Taylor 
series around t = up to 0{t 2 ). As they are only required 
at the values t = ±At and t = ±Ai/2, this should result 
in a scheme converging to second order in At. We obtain 



m a (x,t) = 5 a Al A , 
A(z,*) = 7,a, 
M a A {x,t) = 6 a A + 6 a Bl B A t 



(19) 
(20) 



+ \5 a B (l B cl C A + Ef l)A ) t 2 + 0(t 3 ), (21) 
A A (Z, t) = A (1)A t + ^A {2)A t 2 + 0(t 3 ), (22) 
where we have used the shorthands 



Now we need to address some issues that arise if the 
system of flux-conservative equations under considera- 
tion are a subset of the Einstein equations, in our case 
the Bona-Masso (BM) evolution equations. 

The BM evolution system considers the shift as a given 
function of the space and time coordinates. The lapse can 
be evolved as a dynamical variable, in which case the en- 
tire system is hyperbolic, or it can also be treated as a 
given function of the coordinates, in which case the sys- 
tem is not hyperbolic. In the latter case, one can consider 
the lapse as a given function of the spatial coordinates 
that is constant for one time step, then changes discon- 
tinuously to be a new constant function for the next time 
step. This constant function can then be a functional of 
the Cauchy data at the beginning of the time step. One 
can, for example, obtain the lapse and shift by solving 
elliptic equations for maximal slicing and minimal dis- 
tortion gauge. Such an algorithm has no natural contin- 
uum limit (in time). If one wants to verify convergence 
of the numerical algorithm, one must define an artificial 
continuum limit in which the lapse and shift change at 
intervals At which are multiples of the time step. One 
solves the lapse and shift elliptic equations every time 
step on the coarse grid, every other time step on a grid 
twice as fine, and so on, obtaining a continuum limit in 
which the lapse is constant over finite time intervals. In 
summary, we have three ways of treating the lapse and 
two ways of treating the shift. 

1. Dynamical lapse in the sense of the Bona-Masso 
formalism. There is a hyperbolic evolution equa- 
tion for the lapse and its spatial derivatives. This 
includes "l+log" and harmonic slicing. We shall 
restrict the term "dynamical lapse" to this case. 

2. Non-dynamical lapse or shift that nevertheless de- 
pends on Cauchy data in the manner just described, 
for example maximal slicing. We shall call this a 
"live lapse" or "live shift." 
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3. A non-dynamical lapse or shift that is really a given 
function of space and time coordinates, obtained 
from an analytic solution. This includes the Kerr- 
Schild and related slicings of a single black hole. 
We shall call this an "exact lapse" or "exact shift." 

D. Finite differencing in the x coordinates 

With the differential equation transformed, we now 
have to finite difference in the x grid. This will involve 
creating a new computational grid by interpolation, and 
evolving on that grid rather than the original grid. In 
this section, we explain the details of implementing the 
scheme. 

We excise an irregularly shaped region of the space- 
time from the numerical domain by declaring this set of 
gridpoints to be "masked" . It is technically easier to al- 
locate memory, and to loop over grid points, as if these 
points were still part of the numerical domain, and then 
simply to ignore these points. We do this by overwriting 
these points with, for example, flat spacetime, and set- 
ting a flag. In one dimension, this operation gives no real 
savings in complexity, but in three dimensions, the code 
is significantly simplified by this approach. 

It is useful to first consider the no shift case, where 
causal differencing reduces to ordinary differencing. In 
order to update points neighboring the masked region, 
we have two options. We can evolve all unmasked points 
that depend numerically only on unmasked points, and 
then recreate the remaining unmasked points by extrap- 
olation after the time step. Alternatively, we can evolve 
all unmasked points after first having created any masked 
points they depend on by extrapolation before the time 
step. The work of Scheel et al. Q and the Grand Chal- 
lenge alliance has used extrapolation after the time step. 
We have explored both possibilities, as they can be im- 
plemented with the same code. 

We use cubic interpolation and extrapolation. 

If there is a shift, it will tend to reduce the amount 
of extrapolation needed. Ideally it will turn extrapo- 
lation into interpolation (creating a "boundary without 
boundary conditions"), but in the numerical work dis- 
cussed here, some extrapolation is often needed. In the 
presence of a shift (and therefore a tilt in the stencil) 
all grid points, not just those at the boundary, need to 
be interpolated to new locations. The interpolation and 
extrapolation are dealt with in a single algorithm. To 
demonstrate the location of the grids, in Fig. [| we show 
a one-dimensional Lax-Wendroff stencil with our causal 
differencing approach in the presence of a tilt vector. The 
top picture indicates grid placements when we align grids 
at the beginning of the step, and the bottom indicates 
alignment at the end of the step. 

There is one technical difference between the two al- 
ternatives. When interpolating/extrapolating after the 
time step, we extrapolate only to unmasked points. Be- 



cause the points are unmasked, we know the tilt vector 
at that point, and therefore we know the x coordinate 
values to which we interpolate when re-creating the x 
grid. If we interpolate/extrapolate before the time step, 
we must extrapolate the tilt vector to masked points in 
order to find the location of the x grid in x coordinates 
inside the stencil of the boundary point. 



x x 




x x 



FIG. 2. We show causal Lax-WendrofF stencils, and the 
relationship between the x and x coordinates for the two al- 
ternatives of interpolation/extrapolation before the time step 
(lower diagram) and after the time step (upper diagram). In 
both diagrams, the left stencil shows how a "boundary with- 
out boundary condition" is achieved if the stencil is tilted 
enough; The extrapolation to obtain the leftmost point be- 
comes an interpolation. 

Fig. ^ illustrates the modular construction of the causal 
differencing algorithm. Fig. || assumes interpolation be- 
fore the time step, but an almost identical prescription 
exists for interpolation after the time step. At the begin- 
ning of a time step, all fields are known at point A and all 
other points on the x grid. We interpolate all dynamical 
variables to point C. This is the key step of causal differ- 
encing, namely to interpolate in order to mimic the effect 
of a transport term in the evolution equations. Dynam- 
ical variables are and , plus in the BM formalism 
the Dj k and Vi. Non-dynamical variables are /3 l , and in 
the BM formalism B*. a is non-dynamical in the ADM 
formalism, while in the BM formalism a and Ai can be 
treated either as dynamical or as non-dynamical. 
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t+dt 



t+dt/2 



-0- 



^0- 



-0- 



-dt/2 



-dt 



FIG. 3. Numerical grid showing lines of constant t, con- 
stant x 1 , and constant x l . (Recall that t = t.) Points A and 
E are on the main (x) grid, points B, C and D are not. Sur- 
faces of constant t are labelled on the left by the coordinate 
time, and on the right by the internal time variable used for 
Taylor expansions. At the beginning of the time step, all fields 
are given at point A, at the end of the time step all fields have 
been obtained at point E. 

The non-dynamical variables must be provided at 
points C and D for the Lax-Wendroff algorithm, and at 
points C and E for the MacCormack or MacCormack- 
like algorithms. If these variables are "live" as described 
above, this is done by interpolation. Note that live non- 
dynamical variables are independent of t along lines of 
constant x, so that interpolating these variables to point 
D is the same as interpolating them to B, and they are 
the same at E as at A. If the non-dynamical variables are 
"exact" (derived from an exact solution), they are set at 
C, D and E using the correct value of both x and t. 

Finally, we need to provide M\ and m° for all causal 
schemes, and for the flux-conservative schemes also Aa 
and A, at either C and D (for Lax-Wendroff) or C and 
E (for MacCormack). These are obtained from the tilt 
j A at point E by a Taylor expansion, as described above. 
Recall that by assumption r y A is independent of t along 
lines of constant x, so that j A at point E is the same as 
at point A. In practice Y is a multiple of the shift (3 l . 

If we solve the equations in the x coordinates in the 
flux-conservative form ([ll]), we can use a standard evo- 
lution algorithm, taking into account only the flux and 
source correction terms. Here we have tested a Lax- 
Wendroff and a MacCormack algorithm, both incorpo- 
rating the sources and dealing with all spatial directions 
at once. Note that the presence of flux terms in the source 
correction terms would make Strang splitting these equa- 
tions quite inefficient, and therefore we do not use a 
Strang split. 

Since we will investigate both ( JTl| ) and (|l(]) below, it is 
important to explicitly describe our numerical method for 
the advective form, ([!(]). Differencing the advective form 
with a MacCormack-like method amounts to multiplying 
the finite differences of the fluxes by the appropriate pre- 
factor. That is, terms like -^{F{ui + i) — F(ui)) become 
^M{F(u i+ i) — F(ui)), where M is the appropriate m a 



or M a a factor. 



E. Testbeds 

We shall test our algorithms on the Schwarzschild 
spacetime, in a coordinate system (t, r, 9, <p) adapted to 
spherical symmetry. As the spacetime is static, there are 
many coordinate systems in which all fields are indepen- 
dent of the time coordinate. Without loss of generality, 
we use a radial coordinate r defined so that Awr 2 is the 
area of any surface t = const, r = const. In other words, 
we set gee = r 2 by definition. The remaining coordinate 
freedom is the freedom to slice the spacetime into surfaces 
t = const. One slicing in which all metric coefficients are 
independent of t is of course t = £s c hw> where ischw is 
the usual Schwarzschild time coordinate. This slicing is 
singular at the event /apparent horizon. All other slicings 
which leave the metric coefficients i-independent are of 
the form t — tschw + /(?*)• For one choice of f(r) [which 
diverges at the horizon r = 2M as 2M ln(r — 2M)j one 
obtains the Eddington-Finkelstein, or Kcrr-Schild slicing, 
which is regular at the horizon: 

2M 



a " = g rr = 1 + 

r 

2M ( 2M 

P = ^^7< Kee = 2M[l + 

r + 2M \ r 



-1/2 



2M r + M 

~r 2 ~r + 2M \~ ' r J 



2M\ 
1 + 



1/2 



(26) 



For another choice of /(r), with the same singular part, 
but a different regular part, one obtains the Painleve- 
Gullstrand slicing, in which the 3- metric is flat: 



a = g r 



K r 




Kb 



2M 
r 



V2Mt 



(27) 



Finally, we shall consider the harmonic time-independent 
slicing 



a = g rr = 



4M 2 



( 2AA / 
WM 2 AaM 2 



K r 



AaM 2 ( 3M AM 2 4M 3 \ , . 

2+ + + , (28) 



which has also been used by Scheel et al. Note that 
oT 2 = g rr holds for any time-independent slicing of the 
Schwarzschild spacetime. 
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F. Boundary without boundary conditions 

Fig. |^ illustrates a situation in which we can update 
grid points at the excision boundary without imposing 
any boundary condition and without needing to extrap- 
olate (BWBC). If we align the x and x grids before the 
time step, and the tilt is so large that the entire sten- 
cil lies within the unmasked region, its base points can 
be obtained by interpolation. In the absence of a tilt, we 
would have to extrapolate into the masked region. A sim- 
ilar picture applies when we align the grids after the time 
step, and one sees from the two diagrams that the condi- 
tion on the tilt necessary to obtain a boundary without 
boundary condition is the same for the two alternatives. 

In a spherically symmetric situation, both the tilt and 
the shift are in the radial direction. A simple calculation 
then shows how big the tilt must be to obtain a boundary 
without boundary condition for a spherically symmetric 
black hole. Assuming that both vanish at large radius, it 
is sufficiently general to consider a tilt that is a constant 
multiple of the shift, that is r y A = t/3 a . We call the 
constant r the "tilt factor" . 

We have to take into account two separate conditions: 
The tilt must be large enough to obtain a boundary with- 
out boundary condition, and the physical light cone must 
lie entirely inside the numerical domain of dependence as 
a necessary requirement for numerical stability. Let us 
denote the "Courant number" At/Ar by C, and recall 
that 7 r = t(3 t . Radial null geodesies obey 



(29) 




while the numerical light cones have slopes 



dr 



— = -r(3 r ±C~ 1 . 
at 



(30) 



Therefore the conditions that the inner and outer edge 
of the (past) physical lightcone lie inside the numerical 
lightcone are 



(31) 
(32) 



where 5 is a safety margin (measured in units of Ar). It 
is easy to see that these two conditions are equivalent to 



Cr(3 r - 1 < C yi3 r - 
Crf3 r + 1 > C ( (3 r + 




\(r-lW\< 



1 



c 



(33) 



Note that this condition has to be obeyed on the entire 
grid. The condition that the inner edge of the numerical 
light cone is tilted inwards is easily seen to be 



where e is another dimensionless safety margin. This 
condition needs to be obeyed only at the excision radius, 
and only if we want to avoid extrapolation. 

At large radius, the shift vanishes, while the function 
a l 1 \j9rr determining the width of the light cone is one for 
all usual coordinate systems on Schwarzschild. This gives 
us the stability condition C < 1 — S. In flat spacetime 
this would be all. In Kerr-Schild coordinates, we see that 
a /y/9rr decreases from one via 0.5 at r = 2M to zero 
at r = 0. If we set r = 1 (which was done implicitly in 
previous causal differencing schemes) , the global stability 
condition is simply C < 1 — 8. The shift grows via 0.5 at 
r = 2M to one at r = 0. With r = 1, it never becomes 
quite large enough to allow a BWBC. But a larger value 
of t does allow it, as one can easily see by plotting the 
four slopes (|2f||3(|) against r for different values of C and 

T. 

Conversely, for a given finite differencing method one 
can find the necessary value of t and maximum excision 
radius tq to achieve BWBC. To do this, we take e, S 
and C as given, assume r > 1, set r = r , saturate 
the two inequalities ( |33"| , p4| ) , and solve for r and r . For 
e = 6 = 0.2, C = 0.8 we obtain r = (2/3)M and r = 2. 

In the spatially flat coordinate system, we can gener- 
ally achieve a BWBC quite easily, even and typically with 
r = 1. For e = 5 = 0.2, C = 0.8 we obtain r = (8/9)M 
and r = 1. For 5 — 0.3, C = 0.7 (larger stability margin) 
and e = 0.05 (in one dimension we really need no safety 
margin here) we obtain again = (8/9)M and r = 1, 
and this is the case we have tried numerically. 

In general, the natural choice for S is 6 = 1 — C, which 
makes the stability margin the same at the excision ra- 
dius as at infinity. In this case we obtain t = 2 in 
the Kerr-Schild slicing and r = 1 in the flat slicing of 
Schwarzschild, independently of e. 

In three space dimensions in Cartesian coordinates, the 
shift vector is not typically aligned with a grid axis, and 
therefore factors of up to y/3 arise in various places. It 
is clear that BWBC is then more difficult to achieve. 
Still, it seems possible for the right choice of slicing, a 
tilt factor r > 1, and a sufficiently small excision radius 
ro. Excision of black holes in three dimensions will be 
investigated elsewhere. 



III. THE ONE-DIMENSIONAL BONA-MASSO 
SYSTEM 

A. The equations 

We consider the BM system in spherical symmetry. 
We choose coordinates t and r and a diagonal 3- metric. 
This is the same system as considered in |l2] ], with the 
addition of a shift and conformal factor. 

We have four gauge fields, 



(3 r > 



1 + e 

Ct ' 



(34) 



(a,A r ,/3 r ,B r r ), 



(35) 
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where A r = a, r /a and B r r = (3 r r /2. As discussed 
above, the a and A r can be dynamical, live or exact, 
and the j3 r and B r r are only live or exact. 
We have 7 dynamical variables, 

(ffrr, gee, D rrr , D r ee, K rr , Kee, V r ) , (36) 

where -Dyfc = and V% = -Dy 7 . We note that gee = 

g^ sin 2 9 by spherical symmetry. This has the effect of 
making 



D J 



D, 



D, 



oo 



g r 



and reducing the definition of V r to 

D r ee 



gee 



We also note the useful result 

K rr 



— K J - 



K = K 



gee 



(37) 



(38) 



(39) 



We optionally introduce a conformal rescaling of the 
metric, g — > i/j^g, and define ip r = Vv/V 1 an d VVr = 
ip^r/ifi- The results in exact spacetimes given below will 
not use a conformal rescaling of the metric, but we give 
it for completeness here. 

With these choices the non-vanishing BM sources for 
the Ricci (n = 0) system 113] become 



S-a 
S~g r r 

S-geo 

SJC rr = 2K rr B 



-a 2 fK/iJj 4 +a/3 r A, 
-2aK rr li> A + 4g rr B r r + 

2(3 T D rrr + 4/3 4 -0r5rr, 

-2aK 8 e/i> 4 + 2(3 r D r9e + Wi>r9ee 

K 
g rr 



(40) 

(41) 
(42) 



K rr (K 6 e 



+aA r 



D r 



2ijj r 



2a 



{D r ee + 2ip r gee) ( D 



9r 



SJCes = -2K 0S B, 



gee 

2aA r V r + 8aA r ip r - 
Aa(A r ip r + tf) rr — %p 2 ) 

K rr Kee 



gee 



+ 



Vg r 



(D rrr + 2ip r g rr ){D r ee + 2i\) r gee) 



9rr 



(43) 



1 - 



2a 



A r i> r gee (ip rr - ^ 2 r )ge 



Qrr Qrr 

r ri- 



ff rr 



S.Vr = -2 



[A r K e 



i> 4 gee 
{D r ee + 2i> r gee) 



Kn 



K, 



rr 
Qrr 



(44) 



(45) 



and the non-vanishing fluxes become 

afK/^ - f3 r A r 
aKrr/ip 41 — 2g rr B r r 

(3 D rrr 2/3 tprQrri 

aKee/ip 4 - f3 r D r9 g - 2(3 r ^ r gee, 

D r f)o 



F_A r 



F.D r ee 

F.K rr = -[3 r K rr + a 



F-Kee 
F_V r 



2V r 
D r aa 



A r 



gee 



9r 



-p r (V r + 4Vv 



(46) 

(47) 
(48) 

(49) 

(50) 
(51) 



All other fluxes and sources vanish identically. We note 
that the conformal factor is moved from the fluxes of the 
Kij into the sources, and is not finite differenced, ip and 
its derivatives are given analytically. 
The Ricci scalar of the 3-metric is 



d r D r 



-16 



g rr gee 
D r eeipr 



+ 4- 



D 2 



9rr^ 
^ 'iprr 
Qrr 



J 2 
9rr9oQ 



2 

9ee 



D rrr ip r 

9rr 



g^gee 

and the Hamiltonian constraint, 

Kee{2K rr gee + Kg S g rr ) 



U = R + 2- 



(52) 



(53) 



The maximal slicing equation is 



D r 



Dr 



agr 

V> 4 



Qrr 
K r 
9r> 



2ip r 



Kn 



(54) 



B. Numerical results 

1. E aldington- Finkelstevn 

We present the results of applying our causal differ- 
encing schemes to various black hole spacetimes. As our 
base test, we will use an M — 1 black hole on the domain 
1 < r < 4, so the horizon is at r = 2 and our buffer zone 
has width 1. We use our excision boundary condition at 
the inner boundary, and blend against the analytic solu- 
tion at the outer boundary. We will find, in general, that 
we can keep the Eddington-Finkclstein metric stable for 
many hundreds of M using our schemes. The error in the 
solution grows linearly in time, but converges away faster 
than second order towards zero. In other words, whereas 
Scheel et al. with the unmodified hyperbolic system || 
can achieve run times of order 10M, using the unmodi- 
fied BM system, we can achieve run times in the 100 to 
1000M range. 
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We use the advective and fully flux-conservative causal 
MacCormack algorithms. When we interpolate before 
the step, we use the analytic value of the shift on the in- 
ner boundary point (the only masked point) to re-locate 
the stencil for the first un-masked point. When we in- 
terpolate after the step this is not necessary. We then 
update the first (and only) un-masked point using an ex- 
trapolator, although this is only for a visual effect, and 
does not affect the evolution of the system; with the ex- 
ception of the shift, fields could take any value at the 
innermost (masked) point. 

In all cases we report error as either E = \g rr - 
\gee - 9lT ct \ + \Krr~ K°* act \ + \K ee - K%% act | or as a norm 
over Hamiltonian constraint violation. Both properties 
show the same convergence behavior. The norm 1 1 is the 
sum of absolute values divided by the number of grid 
points. 

In Fig. U we see the Hamiltonian constraint for the 
Finkclstein black hole evolved with the advective scheme. 
The error grows smoothly until it is of order one, and the 
code crashes. Fig. ^ below shows that the growth of the 
error is between linear and quadratic for the first 200AZ 
(at the highest resolution). Later, the growth leading to 
the crash is approximately exponential. 



^cxact I 
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FIG. 4. We show the evolution of the norm of the Hamilto- 
nian constraint for a single Eddington-Finkelstein sliced black 
hole evolved with our advective MacCormack-like scheme. We 
note that generically runtimes are long and errors are small 
until the code crashes. The convergence rate implied by this 
graph is about 2.7 until the system crashes. Note that M = 1 
throughout this paper. 

We repeat this comparison in Fig. ^ using the fully 
flux-conservative MacCormack scheme, and see the same 
behavior. That is, we see that the error converges to zero 
faster than second order, and runs are generally cut short 
by a crash as the error grows faster-than-linearly towards 
the end of the run. 



0.8 




1000 1200 



FIG. 5. We show the Hamiltonian constraint violation at 
three different resolutions for the Eddington-Finkelstein black 
hole. Compare with Fig. ^j. 

Despite this initial transient, the Hamiltonian con- 
straint violation at late times are essentially the same in 
both schemes. In Fig. ^ we compare the error at a fixed 
resolution between the two methods, and see exactly this 
behavior. 



0.020 



0.015 



0.000 



Non-flux conservative 




200 



FIG. 6. We compare the error at Ar = 3/200 for the fully 
flux conservative and advective methods. We note that, after 
an initial transient error in the fully flux-conservative scheme, 
the errors are essentially the same, and both errors grow in 
time. 

The source of the error in our solution is almost en- 
tirely dissipation of the solution at the inner boundary. 
In Fig. [?] we show plots of g rr on the Ar = 3/200 grid at 
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times OA/, lOOAf, and 200Af . We note that the entire er- 
ror comes from dissipation at the inner boundary, where 
the solution slowly "melts" away. Our scheme does not 
iterate towards a static solution, it seems, but rather g rr 
drops linearly (but very slowly) in time, as shown in the 
inset. 



3.0 




FIG. 7. We show the evolution of g rr for BM with advective 
MacCormack using 200 points (Ar = 3/200). In the large 
figure we show g rr at 0, lOOAf, and 200A4". In the inset we 
show the dissipation of g rr at the inner boundary, which is 
the princiapl source of error. 



o.io r 



o.oo 



1 

S -0.10 

1 

3 

a 



-0.20 



0.0020 
0.0015 



g 0.0010 
o 



g 0.0005 

a 

0.0000 




-0.0005 L_ 
2.0 



2.5 



3.0 



3.5 



4.0 



-0.30 t 



FIG. 8. We show the Hamiltonian constraint for the BM 
system with the advective causal MacCormack using 200 
points (Ar = 3/200) shown at t = 0, t = lOOAf, and 
t — 200Af. We notice that the Hamiltonian constraint vi- 
olation (which converges to zero at second order) is large 
inside the horizon (r < 2M), but in the area outside the 
horizon (r > 2M) shown in the inset, the violation has not 
"escaped." This figure confirms one of the fundamental ideas 
of the AHBC, namely, that an inaccurate but stable scheme 
applied inside the horizon will not affect the system outside 
the horizon. 



One of the promises of the AHBC is that errors made 
at the boundary inside the horizon won't affect evolu- 
tion outside the horizon. In Fig. |^ we show this effect. 
We measure the Hamiltonian constraint at a given reso- 
lution (again, Ar = 3/200). We observe the constraint 
converging towards zero at or above second order, and 
therefore can be concerned with its magnitude. Notably, 
we can see that the violation is several orders of magni- 
tude larger inside the horizon. This gives us hope that 
our technique is correctly obeying the causal structure of 
our spacetime; numerical errors in the exterior are barely 
affected by (large but stable) dissipative boundary errors 
in the interior. 



2. Spatially flat Schwarzs child, and boundary without a 
boundary condition 

We have tested all our numerical schemes on a second 
slicing of Schwarzschild, the spatially flat one discussed 
above. For a direct comparison with the Eddington- 
Finkelstein slicing, we have used the same Courant factor 
C = 0.5 and excision radius ro = M. Results are very 
similar. Our main observations here are that the finite 
differencing schemes again do not tend to a static so- 
lution of the finite difference equations, the error grows 
approximately linearly in time, and a doubling of spatial 
resolution therefore buys a run time that is four times 
as long. Comparison with the Eddington-Finkelstein re- 
sults confirms the point that our algorithm is not spe- 
cially made for a particular solution, works better with 
higher resolution, and is therefore generic and robust. 
Again, results are much better, at the same resolution, 
with causal differencing than without. 

In a second series of runs, we have tested a param- 
eter choice, C = 0.7 and ro = 0.9Af, that allows us 
to obtain an excision boundary without extrapolation. 
To our knowledge this is the first time that black hole 
excision was achieved using causal differencing without 



10 



extrapolation at the excision boundary. Nevertheless, 
the experimental result is unspectacular: only somewhat 
larger runtimes are achieved. This negative result is im- 
portant as it indicates that extrapolation at the excision 
boundary is not the prime cause of error and eventual 
instability. Having a no-extrapolation boundary may be- 
come more useful in three space dimensions, where three- 
dimensional extrapolation on a jagged "Lego" boundary 
is not as straightforward as in spherical symmetry. 

3. Harmonic slicing 

As pointed out by Scheel et al. 0, there is a time- 
independent slicing of Schwarzschild that is also har- 
monic. In one set of runs, we have used this as our third 
slicing using an exact lapse and shift. In a second set 
of runs, we have evolved the same initial data with a 
live lapse, namely the harmonic slicing condition. With 
the exact lapse and shift, we find again the same scaling 
and run times as for the other two slicings. With the 
live harmonic slicing lapse, runtimes at high resolutions 
are in excess of 22000 M, much longer than for the ex- 
act lapse: here the deviation from the true solution does 
not increase monotonously, but turns around and settles 
down. As this turnaround occurs at large deviations, 
these longer run times are essentially accidental. Con- 
vergence at small errors, however, is still second order, 
and the live lapse seems to be as stable as the exact one 
we used in the other runs. 

In order to obtain a more direct comparison with 
the results of Scheel et. al 0, we have made runs 
with exactly their resolution, setting the excision radius 
much closer to the horizon (at r = 1.8M, instead of 
l.OJVf) and/or the outer boundary much further out (at 
r ~ 120M instead of 4.0M). The combination inner 
boundary close to the horizon - outer boundary close in 
is numerically unstable. (We do not know why.) The 
combination inner boundary close to the horizon - outer 
boundary far out is again stable. The combination in- 
ner boundary at 1.0M - outer boundary far out works 
equally well. In summary, in similar circumstances we 
obtain similar results as Scheel et al., although the two 
codes use different hyperbolic formulations of the Ein- 
stein equations. 

Our runs are summarized in Table 1. 



In no case does the solution of the finite difference equa- 
tions settle down to stable state, so that all runs crash 
after a finite time. Nevertheless, the numerical error is 
as well-behaved as one can hope for: it is proportional to 
h 2 on the one hand, where h is the numerical resolution 
in space and time, and while it is small, it grows between 
linearly and quadratically with time on the other hand. 
This means, and experience confirms, that with twice 
the number of radial grid points one roughly doubles to 
quadruples the run time before crashing. 

Our results are similar to those of Scheel et al. 
In both cases, causal differencing was applied to a hyper- 
bolic formulation of the Einstein equations. Both inves- 
tigations find the same behavior of the numerical error. 
At the same resolution, Scheel et al. in their more recent 
work |0 report run times about a factor of 10 larger than 
ours. One of our four causal differencing schemes (advec- 
tive with interpolation at the end) is similar to that of 
Scheel et al. The differences are as follows. We have used 
an exact shift, and either an exact lapse or the harmonic 
slicing condition, while Scheel et al. use the harmonic 
slicing condition and live minimal distortion shift. We 
excise at a fixed coordinate radius, while Scheel et al. at- 
tach their excision radius to the apparent horizon as it 
changes through numerical error. In spherical symmetry, 
these differences are probably not as important as the 
finite differencing scheme itself. 

On physical grounds, no boundary condition is re- 
quired at the excision boundary - it is not a timelike 
boundary, but a future spacelike one. (The term "ap- 
parent horizon boundary condition" is misleading in this 
sense, and one better speaks of "black hole excision".) 
We have used causal differencing to obtain a genuine 
boundary without boundary condition, that is, without 
numerical extrapolation at the excision boundary. This 
works well, but does not seem to have a numerical ad- 
vantage, at least in spherical symmetry. In other words, 
the extrapolation boundary condition does not appear 
to be the dominant cause of numerical error in spheri- 
cal symmetry. (This may well be different in three space 
dimensions.) 

Our causal differencing methods are immediately ap- 
plicable to the Bona-Masso formulation of the Einstein 
equations in three dimensions. Ongoing work on black 
hole excision in three dimensions will be reported else- 
where. 



IV. CONCLUSIONS 

Our main result is that while non-causal differencing, 
with an extrapolation boundary condition, does not work 
for black hole excision, all causal schemes do. The idea 
of causal differencing is robust in the sense that all four 
schemes we have implemented perform similarly well, and 
that no modifications of the field equations were required. 
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Method 








Runtime in units ol M 


Causal 


Form 


At/ A 


r Interpolate 




Low 


Med 


High (Ano„ /A) 


Eddington-Finkelstein slicing, Ari ow = 0.06, ro = M, r max = AM 


N 


FC 


0.5 






28 


47 


95 


N 


FC 


0.2 






10 


16 


45 


N 


FC 


0.1 






7 


8 


12 


C 


Adv 


0.5 


start 




118 


390 


1412 


C 


Adv 


0.5 


end 




66 


170 


429 


c 


FC 


0.5 


start 




100 


339 


1249 


c 


FC 


0.5 


end 




58 


152 


389 








Spatially flat slicing, Ari ow 


= 0.06, r = M, r max = AM 




N 


FC 


0.5 






8.8 


7.5 


1.2 


c 


Adv 


0.5 


start 




76 


293 


1165 


c 


Adv 


0.5 


end 




45 


146 


451 


c 


FC 


0.5 


start 




49 


293 


787 


c 


FC 


0.5 


end 




31 


107 


339 






Spatially flat slicing, Ari ow = 


0.06, r 


= 0.9M, r max = AM 




N 


FC 


0.7 






0.6 


0.3 


0.2 


c 


Adv 


0.7 


start 




127 


324 


1411 


c 


Adv 


0.7 


end 




35 


82 


503 


c 


FC 


0.7 


start 




33 


26 


42 


c 


FC 


0.7 


end 




21 


32 


33 








Harmonic slicing, An ow = 


0.06, r 


= M, r max = AM 




N 


FC 


0.5 






16 


29 


38 


N 


FC 
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11 
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c 
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start 
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247 
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c 
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end 
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c 


FC 
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start 
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c 


FC 


0.5 


end 
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Live harmonic lapse, harmonic slicin 


5, An ow 


= 0.06, r = M, r max = AM 




N 


FC 


0.5 






25 


51 


109 


N 


FC 


0.1 






15 


20 


31 


c 


Adv 


0.5 


start 




175 


> 22000 


> 22000 


c 


Adv 


0.5 


end 




121 


> 22000 


> 22000 


c 


FC 


0.5 


start 




156 


497 


> 22000 


c 


FC 


0.5 


end 




116 


3370 


> 22000 






Live harmonic lapse, harmonic slicing, 


An ow = 


0.125, r = 1.8M, r max = AM 




c 


Adv 


0.5 


start 




64 


68 


68 


c 


Adv 


0.1 


start 




53 


55 


55 






Live harmonic lapse, harmonic 


dicing, 


Ariow = 


0.125, r = M, r max = 120M 




c 


Adv 


0.5 


start 




16 




2103 


Live harmonic lapse, harmonic slicing, An ow = 0.125, r = 1.8M, r max = 120M 


C 


Adv 


0.5 


start 




68 


871 


1389 



TABLE I. Summary of black hole evolutions in one dimension. Causal is either non-causal with extrapolation boundary 
condition (N) or causal (C). Form of the equations is either flux-conservative (FC) or advective (Adv). The integration method 
is MacCormack or MacCormack-like (for the advective form). Interpolation is either not needed (-), or done at the start or 
the end of the evolution step. Note that causal advective MacCormack with interpolation after the time step is essentially the 
method of Scheel et al. Note that plain (non-causal) MacCormack, with the same extrapolation boundary that works for the 
causal schemes, is unstable in all situations tried, even for a very small Courant number. 
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